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ABSTRACT 

We describe techniques for incorporating feedback from star formation and black hole 
accretion into simulations of isolated and merging galaxies. At present, the details of 
these processes cannot be resolved in simulations on galactic scales. Our basic approach 
therefore involves forming coarse-grained representations of the properties of the in- 
terstellar medium and black hole accretion starting from basic physical assumptions, 
so that the impact of these effects can be included on resolved scales. We illustrate 
our method using a multiphase description of star-forming gas. Feedback from star 
formation pressurises highly overdense gas, altering its effective equation of state. We 
show that this allows the construction of stable galaxy models with much larger gas 
fractions than possible in earlier numerical work. We extend the model by including a 
treatment of gas accretion onto central supermassive black holes in galaxies. Assuming 
thermal coupling of a small fraction of the bolometric luminosity of accreting black 
holes to the surrounding gas, we show how this feedback regulates the growth of black 
holes. In gas-rich mergers of galaxies, we observe a complex interplay between star- 
bursts and central AGN activity when the tidal interaction triggers intense nuclear 
inflows of gas. Once an accreting supermassive black hole has grown to a critical size, 
feedback terminates its further growth, and expels gas from the central region in a 
powerful quasar-driven wind. Our simulation methodology is therefore able to address 
the coupled processes of gas dynamics, star formation, and black hole accretion during 
the formation of galaxies. 

Key words: galaxies: structure - galaxies: interactions - galaxies: active - galaxies: 
starburst - methods: numerical 



1 INTRODUCTION 

It is now recognised that galaxy collisions and mergers are 
relevant to a wide range of phenomena associated with both 
ordinary and active galaxies. The collisionless dynamics of 
this process is relatively well-understood. The seminal work 
of Toomre & Toomre (1972) using restricted N-body meth- 
ods demonstrated that gravitational tidal forces between 
disks interacting transiently can account for the narrow 
tails and bridges commonly seen in morphologically pecu- 
liar galaxies. Subsequent work using self-consistent models 
(Barnes, 1988, 1992; Hernquist, 1992, 1993c) verified these 
conclusions and made it possible to extend the calculations 
to investigate mergers and the detailed structure of the rem- 
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nants left behind. These studies showed that in many re- 
spects the remnants closely resemble elliptical galaxies, as 
proposed by Toomre (1977). 

However, it is clear that pure stellar dynamics is not 
adequate to explain the extreme variety of objects linked to 
galaxy encounters. This possibility was already anticipated 
by Toomre & Toomre (1972) who suggested that collisions 
could "bring deep into a galaxy a fairly sudden supply of 
fresh fuel," leading to a period of violent, enhanced star 
formation. The significance of non-gravitational processes 
in galaxy interactions was confirmed in the 1980s when it 
was recognised that the brightest infrared sources seen with 
the IRAS satellite are invariably associated with peculiar 
galaxies (Sanders et al., 1988; Melnick & Mirabel, 1990). A 
natural interpretation of these systems is that the infrared 
emission is powered by a starburst, triggered through merg- 
ers. 

There are also numerous studies indicating that 
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quasars, radio galaxies, and active galactic nuclei (AGN) 
are preferentially found in tidally disturbed objects (for re- 
views, see e.g. Barnes & Hernquist, 1992; Jogee, 2004). In- 
deed, much circumstantial evidence suggests an evolutionary 
connection between merger-induced starbursts and quasar 
activity (Sanders et al., 1988). In addition, in recent years, 
strong dynamical evidence has been accumulated indicating 
that supermassive black holes reside at the centre of most 
galaxies (Kormendy & Richstone, 1995; Magorrian et al., 
1998; Kormendy & Gebhardt, 2001). Moreover, a remark- 
able connection between supermassive black holes and the 
properties of their host galaxies has been discovered: the 
masses of central supermassive black holes arc tightly corre- 
lated with the stellar velocity dispersion of their host galaxy 
bulges (Fcrrarcsc & Mcrritt, 2000; Gebhardt et al., 2000; 
Tremaine et al., 2002), as well as with the mass of the 
spheroidal component (Magorrian et al., 1998; McLure & 
Dunlop, 2002; Marconi & Hunt, 2003). This strong hnk sug- 
gests a fundamental connection between the growth of black 
holes and the formation of stellar spheroids in galaxy halos 
(Kauffmann & Haehnelt, 2000; Volonteri et al., 2003; Wyithe 
& Loeb, 2003; Granato et al., 2004; Haiman et al., 2004; Di 
Matteo et al., 2003, 2004a, and references therein). 

A better understanding of galaxy interactions and 
mergers will, therefore, require simultaneous accounting of 
the physics of interstellar gas, star formation, the growth 
of supermassive black holes in galactic nuclei, and various 
forms of feedback associated with both massive stars and 
AGN. As we discuss in detail below, it is not immediately 
obvious how these effects should be incorporated into sim- 
ulations, primarily because we presently lack a sufficiently 
developed theory of star formation, but also because the 
current generation of computer models cannot resolve the 
complex structure of star-forming gas on the scales of whole 
galaxies. For these reasons, efforts to study non-gravitational 
processes in galaxy mergers have to rely on strong simplifi- 
cations, which have often been very restrictive in previous 
work. 

Negroponte & White (1983) used a sticky-particle treat- 
ment to show that galaxy mergers can concentrate gas in 
the inner regions of a remnant. Later, Hernquist (1989) and 
Barnes & Hernquist (1991, 1996) examined the fate of gas in 
minor and major mergers with a smoothed particle hydro- 
dynamics (SPH) algorithm, neglecting star formation and 
feedback and approximating the structure of the interstellar 
medium (ISM) with an isothermal equation of state. The 
first serious attempts to model star formation in galgixy 
mergers were made by Hernquist & Mihos (1995) and Mihos 
& Hernquist (1996), who also used an SPH method with 
an isothermal gas and included a minimal form of kinetic 
feedback from massive stars (Mihos & Hernquist, 1994b). 
Gerritsen & Icke (1997) and Springel (2000) began to gener- 
alise these calculations by implementing simple approaches 
to allow for departures from an isothermal equation of state. 

While the earlier calculations yielded many successes by 
providing an explicit theoretical link between galaxy merg- 
ers and starbursts, many related phenomena remain poorly 
understood, partly because of the approximations made in 
handling the consequences of feedback. In this paper, we in- 
troduce a new methodology to incorporate star formation 
and black hole growth into galaxy-scale simulations using 
a sub-resolution approach. Starting from a "microscopic" 



physical theory for e.g. the ISM, we form a "macroscopic" 
coarse-grained representation that captures the impact of 
star formation and black hole growth on resolved scales. This 
strategy is not tied to a particular underlying model and can 
be used to graft any well-specified theoretical description 
onto the simulations. As a specific example, in what follows 
we will illustrate our technique using a sub-resolution mul- 
tiphase model for star-forming gas developed by Springel 
& Hernquist (2003a). As an additional component, we will 
add a treatment of gas accretion onto supermassive black 
holes, modelled as coUisionless 'sink' particles, and we treat 
feedback processes associated with the accretion. 

Note that the use of a sub-resolution approach to in- 
clude unresolved effects in simulations is a general concept 
also used, for example, in N-body representations of coUi- 
sionless fluids, where six-dimensional phase space is discre- 
tised into fluid elements that are computationally realised 
as particles. Our treatment is analogous to this, except that 
we endow the particles with internal degrees of freedom that 
characterise, e.g., the thermodynamic state of star- forming 
gas. These internal degrees of freedom evolve in a manner 
consistent with the basic theory from which they are de- 
rived, enabling us to couple the "microphysics" of the ISM 
to the larger scale dynamics of the galaxies. 

This paper is organised as follows. We first discuss the 
construction of compound galaxy models used in our study 
in Section 2. We then summarise our microscopic descrip- 
tions of star formation and black hole growth in Sections 3 
and 4, respectively, and discuss aspects of our numerical 
treatment in Section 5. We use our models to examine the 
impact of various forms of feedback on isolated disks in Sec- 
tion 6 and major mergers in Section 7. Finally, we conclude 
in Section 8. 



2 MODEL GALAXIES 

Each galaxy in our study consists of a dark matter halo, a 
rotationally supported disk of gas and stars, and a central 
(optional) bulge. The parameters describing each compo- 
nent are independent, so that a wide range of morpholog- 
ical types can be specified. The models are constructed in 
a manner similar to the approach described by Hernquist 
(1993a) and Springel (2000), but with a number of refine- 
ments. In particular, we use halo density profiles motivated 
by cosmological simulations, and include pressure gradients 
in setting the equilibrium structure of the gas. The former 
consideration is important for the development of tidal fea- 
tures in interacting systems (Dubinski et al., 1996, 1999; 
Mihos et al., 1998; Springel & White, 1999), while the latter 
improvement is significant when the gas fraction of the disk 
is large, or if the gas has a relatively stiff equation of state. 

2.1 Density profiles 
2.1.1 Halos 

In the earlier work of Hernquist (1993a), halos were de- 
scribed by truncated isothermal spheres with constant den- 
sity cores. However, N-body simulations such as those of Du- 
binski & Carlberg (1991) and Navarro et al. (1996) (here- 
after NFW) indicate that cosmological halos have density 
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Figure 1. Density profiles of NFW and Hernquist model halos, 
matched to each other as described in the text. The halo has 
concentration c = 10. The total mass of the Hernquist model is 
equal to he mass of the NFW profile within the virial radius r2oo ■ 



profiles that are centrally peaked and drop at large radii 
more rapidly with radius than isothermal spheres. 

In view of this, we model the dark matter mass distri- 
bution with a Hernquist (1990) profile, i.e. 

, , Mdm a 



Pdm[ 



2n 



(1) 



with cumulative mass profile M{< r) = Mdm''^/(»' + o,)^- 
There are several motivations for this choice. In its inner 
parts, the shape of this profile agrees with the NFW fitting 
formula, but due to its faster decline in the outer parts, the 
total mass converges, allowing the construction of isolated 
halos without the need for an ad-hoc truncation. Further- 
more, in many situations it is convenient to work with com- 
ponents having analytical distribution functions, as is the 
case for the Hernquist (1990) profile, but not for the NFW 
model. 

To make contact with common descriptions of halos in 
cosmological simulations, we associate the Hernquist profile 
with a corresponding NFW-halo with the same dark matter 
mass within r2oo, the radius at which the mean enclosed 
dark matter density is 200 times the critical density. We also 
require that the inner density profiles are equal (i.e. pdm = 
Pnfw for r <^ ^200), which implies a relation between a and 
the scale length of the NFW profile. The latter is often 
given in terms of a concentration index c, conventionally 
defined as c = r2oo/r3, where is the scale length of the 
NFW halo. We then have the relation 



■rsy/2[ln{l + c)~c/{l + c)]. 



(2) 



For a concentration of c = 10, this gives a ~ 1.73 rs. The 
Hernquist profile then contains 75% of its total mass within 
r2oo, and 99% within 1.1 r 200- 

In Figure 1, we compare the density profile of an NFW 
halo with a Hernquist model of the same concentration. In 
the inner regions, the two mass distributions match closely, 
while at large radii the density profiles asymptote to Pnfw oc 

and /OHcm oc r~^, respectively. At present, observations 



indicate that either form can give a good description of the 
density profiles of actual halos when extended beyond the 
virial radius (Rines et al., 2000, 2002, 2003, 2004). 



2.1.2 Disks and bulges 

We model disk components of gas and stars with an expo- 
nential surface density profile of scale length h; i.e. 



Sgas(»') " 

S*(r) = 



exp(-r//i), 



exp{-r/h), 



(3) 



(4) 



so that the total mass in the disk is Md = (Mgas -|- M*) = 
TidMtot, where nid is dimensionless and Mtot is the total 
mass of the galaxy. 

The bulge is taken to be spherical, for simplicity. We 
also model it with a Hernquist profile: 

(5) 



Pb(0 = -TT- 



2tv r{r + b)3' 

We treat the bulge scale length 6 as a free parameter that 
we parameterise in units of the disk scale length, and specify 
the bulge mass as a fraction rrib of the total mass, i.e. Mb — 
m.bA'/tot. 

We set the disk scale length h by relating it to the an- 
gular momentum of the disk. Following Mo et al. (1998) 
and Springel & White (1999), we first write the total halo 
angular momentum of the halo as 



J 



\ ^1/2 Af3/2 1/2 / 2 



1/2 



(6) 



where A is the usual spin parameter, which we take as a free 
parameter of our galaxy models. The factor 



fc 



:[l-l/(l + c)^-2 1n(l + c)/(l + c)] 
2[ln(l + c) -c/(l + c)]'' 



(7) 



depends only on the halo concentration index. Of direct 
relevance for the structure of the galaxy is the fraction of 
the angular momentum in the disk. We typically assume 
Jd = mdJ, corresponding to conservation of specific angu- 
lar momentum of the material that forms the disk. Assuming 
that the disk is centrifugally supported, this then implies a 
one-to-one relation between A and h. In order to obtain this 
connection, we compute the disk angular momentum assum- 
ing strict centrifugal support of the disk and negligible disk 
thickness compared to its scale-length. Then we can write 



Jd 



Md/ K(i?)(f) exp(_f) 



d7?, 



where the circular velocity is given by 
G[Md„,(<i?)-hMb(<i?)] 



R 

X [Io{y)Ko{y) - 



+ 



2GAfd 2 

— I — y 



(8) 



(9) 



Here y = R/{2h), and the /„ and A"„ are Bessel functions. 
Note that we will later drop the thin-disk approximation for 
the disk's potential in favour of an accurate representation 
of a thick disk. However, for the conversion of the free pa- 
rameter A into a disk scale length h, the accuracy of equation 
(9) is sufficient. 
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We specify the vertical mass distribution of the stars in 
the disk by giving it the profile of an isothermal sheet with 
radially constant scale length zo. The 3D stellar density in 
the disk is hence given by 

, 9 / 2 \ / -R^ 
'h. 



■ sech 



(10) 



We treat zq as a free parameter that effectively deter- 
mines the 'temperature' of the disk, and set the velocity 
distribution of the stars such that this scale-height is self- 
consistently maintained in the full 3D potential of the galaxy 
model. 

However, a similar freedom is not available for the 
gas, because once cooling and star formation processes are 
taken into account, we cannot choose the temperature freely. 
Rather, in a broad class of models the gas will stay close to 
an (effective) equation of state of the form P — P{p), imply- 
ing a tight relation of gas pressure and gas density. Examples 
for such situations include isothermal gas, or the subresolu- 
tion multiphase model by Springel & Hernquist (2003a). For 
a given surface density, the vertical structure of the gas disk 
then arises as a result of self-gravity and the pressure given 
by this equation of state, leaving no freedom for prescrib- 
ing a certain vertical profile. In this situation, the vertical 
structure has to be computed self-consistently, as discussed 
next. 



2.1.3 Vertical structure of the gas disk 

We assume that the vertical structure of the gas disk in 
our axisymmetric galaxy models is governed by hydrostatic 
equilibrium, i.e. 



Pg dz dz 



(11) 



where "1> is the total gravitational potential of all mass com- 
ponents. This equation can be rewritten as 



dz 



7P dz ' 



(12) 



where 7 = is the local polytropic index of the equation 

of state. Note that for a given potential $, the solution of 
this equation is determined by the integral constraint 



^gas(-R; z) 



Pg{R,z) dz, 



(13) 



where Egaa(r) is the surface mass density we prescribed in 
equation (3). Assuming a given potential for the moment, 
we solve equations (12) and (13) by integrating the differen- 
tial equation for a chosen central density value in the z = 
midplane. This chosen starting value is then adjusted until 
we recover the desired surface density for the integrated ver- 
tical structure of the gas layer. We carry out this process as 
a function of radius, using a fine logarithmic grid of points 
in the R— z plane to represent the axisymmetric gas density 
distribution. 



2.2 Evaluating the potential 

The above invokes the problem of determining the potential 
and the resulting gas distribution self-consistently. We ad- 
dress this problem by an iterative method. Starting with 
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Figure 2. Time evolution of a Hernquist dark matter halo with 
initial conditions constructed in the way described in this paper. 
The initially realised input profile is the thin red line. Because 
the approximation of a locally Gaussian velocity distribution is 
not exact, the central profile is not in perfect equilibrium in the 
beginning. As a result, the density in the centre fluctuates down- 
wards on a short timescale of order the crossing time, and then 
relaxes to a slightly softer central profile. After time t = 0.08 
(in units of the dynamical time, idyn = ''200 /f 200, of the halo), 
the deviation is close to its maximum (dashed line). Already at 
t = 0.16, however, a stable profile is reached, which then remains 
essentially invariant with time, as illustrated by the further out- 
put times that are overplotted (at times t = 0.24, 0.32, 0.48, and 
0.64). The vertical dotted line marks the scale below which the 
force law is softened compared to the Newtonian value. 



a guess for an initial potential, we compute the implied 
vertical gas structure as described above, treating the po- 
tential as fixed. Then, we recompute the potential for the 
updated gas distribution, and repeat the whole process. In 
doing this, the potential and the gas distribution converge 
quickly, within a few iterations, to a self-consistent solution. 

However, in doing this, we also need a method to accu- 
rately compute the potential and the resulting force field for 
the radially varying density stratification of the gas. Since 
we want to be able to set up quite massive gas disks that are 
in equilibrium right from the start, we cannot utilise simple 
thin-disk approximations, but rather need to compute the 
full potential accurately. The generality of the gas distribu- 
tions used here makes this a non-trivial exercise. 

We use a somewhat contrived, yet extremely flexible so- 
lution for this task. We set up a discretised mass distribution 
that represents the desired disk components accurately, and 
then compute the potential numerically with a hierarchical 
multipole expansion based on a tree code. For a given (cur- 
rent) mass model of the galaxy, we use for this purpose a 
suitably distorted grid of particles of typically 2048 x 64 x 64 
particles per component (the different numbers refer to ra- 
dial, azimuthal, and vertical directions, respectively). While 
we find that this gives sufficient accuracy for our purposes, 
we note that the number of particles can be easily made 
much larger to generate a still smoother potential, if needed. 
Unlike in a normal N-body code, we are here not interested 
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Figure 3. Rotation curves of model galaxies with the following parameters: (a) 11200 = 160 km s~^, c = 9, = 0.041, mi-, = 0, 
/gas = 0.1, A = 0.033, Jd = 0.041, ZQ = 0.2; (b) same as model (a) but with mb = 0.01367. The resulting disk scale lengths are 
2.74/i^^kpc and 2.46/i~^kpc, respectively. 



in evaluating the forces or the values of the potential for all 
these particles. We only use them as markers for the mass, 
and instead compute the potential at a fixed set of spatial 
coordinates. The work required for this potential compu- 
tation then scales only with log(A''). The primary cost of 
making A'^ large lies therefore in the memory consumed on 
the computer, and not in the CPU time. The potential and 
force field of the spherical mass distributions used for the 
dark matter halo and the bulge are known analytically, so 
we simply add them to the result obtained with the tree for 
the disk components. 



2.3 Velocity structure 

Once the full density distribution has been determined, we 
compute approximations for the velocity structure of the col- 
lisionless and gaseous components. For the dark matter and 
stars in the bulge, we assume that the velocity distribution 
function only depends on energy E and the Lz-component 
of the angular momentum. Then, mixed second order mo- 
ments of the velocity distribution vanish, (vrVz) = (vzV^) = 
{vrv^} = 0, as well as the first moments in radial and ver- 
tical directions, {vr) = (vz) = 0. The velocity distribution 
can then be approximated as a triaxial Gaussian with axes 
aligned with the axisymmetric coordinate system. The non- 
vanishing second-moments can be obtained with the Jeans 
equations. We have 



(14) 



where p is the density of the mass component under consid- 
eration. For the azimuthal direction we have 



/ 2\ / 2\ ^ R d{p{v%)) , 2 

{^d,} = [Vr) H ^„ + Vc 



dR 



where 
2 



(15) 



(16) 



is the circular velocity. 

In the azimuthal direction, there can be a mean stream- 
ing component (v^), which is not determined by the Jeans 
equations. Once it is specified, the dispersion of the Gaus- 
sian velocity distribution in the azimuthal direction is given 
by 



(4) - {'"'t'f 



(17) 



For the stellar bulge, we set (v^) to zero, meaning that the 
bulge is assumed to have no net rotation. For the dark mat- 
ter halo, we assume that it has the same specific angular 
momentum as the disk material. For simplicity, we impart 
this angular momentum by making (^0) a fixed fraction of 
the circular velocity, i.e. {«</,) — fsVc (Springel & White, 
1999). Note that we then have fs <^ 1, i.e. the net stream- 
ing of the resulting dark halo is quite small. 

The velocity structure of the stellar disk can in principle 
be much more complicated, and it will in general have a 
distribution function that does not only depend on E and 
Lz. For simplicity, we however continue to approximate the 
velocity distribution with a triaxial Gaussian that is aligned 
with the axisymmetric coordinated vectors. 

As Hernquist (1993a) points out, observationally there 
is good evidence that (ujj) is proportional to (w?) for the 
stellar disk. We hence assume aR = {vii) = fR{y1^. Note 
that larger values of /_b will increase the Toomre Q parame- 
ter, making the disk more stable against axisymmetric per- 
turbations. In most of our default models, we actually set 
/fl = 1 which corresponds to the approximations made for 
the dark halo and the bulge. Note that the Toomre Q also 
varies in response to the vertical dispersion (w?) which sen- 
sitively depends on the assumed disk thickness. 

To specify the mean streaming, we employ the epicyclic 
approximation, which relates the radial and vertical disper- 
sions in the stellar disk: 



„2 

2 _ dR 



(18) 
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where 

~ RdR \RdR dRy ■ ^ ' 

Using equation (15) for (w^), we then set the streaming 
velocity for the stellar disk as 

M={{-i)-$y'- (20) 

This completes the specification of the velocity structure of 
the coUisionless components. 

For the gas, we deal with a single valued velocity field, 
where only the azimuthal streaming velocity has to be spec- 
ified. The latter is dotorminod by the radial balance between 
gravity on one hand, and centrifugal and pressure support 
on the other hand. Specifically, we here have 

Note that the gas temperature at each coordinate is given 
by the equation of state based on the self-consistent mass 
distribution derived above. 

We carry out the integrations over the force field needed 
in the Jeans equations with the help of a fine logarithmic 
grid in the _R-z-plane. This grid is also used to tabulate the 
final velocity dispersion fields. The difi'ercntiations needed 
in equation (15), for example, are approximated by finite- 
differencing off this grid. Once all of the density and veloc- 
ity distributions have been computed, we initialise particle 
coordinates and velocities by drawing randomly from the 
respective distributions. Values for the velocity structure at 
individual particle coordinates are obtained by bilinear in- 
terpolation of the R-z grid. 

Kazantzidis et al. (2004) recently examined the accu- 
racy of initialising spherically symmetric dark matter halos 
with a Gaussian velocity distribution with dispersion derived 
from the Jeans equations. They found that the innermost 
parts of halos are not precisely in equilibrium when this ap- 
proximation is applied. As a result, the core relaxes within 
the first few crossing times to a density profile which lies 
slightly below the input profile in the centre, such that the 
central cusp becomes softer. Using a more accurate compu- 
tation of the velocity distribution function that takes higher 
order moments into account, they also demonstrated that 
this behaviour can in principle be avoided. 

In Fig. 2, we show an example of this effect. To illus- 
trate this in a manner unaffected by two-body relaxation we 
have set-up a pure dark matter halo with 4 million particles 
using our code for constructing compound galaxies. When 
evolved forward in time, the density profile in the very cen- 
tre fluctuates below the desired input value during the first 
crossing times. Afterwards, the density profile shows no fur- 
ther secular evolution and stays extremely stable for times 
of order the Hubble time. The magnitude of the initial relax- 
ation effect appears to be consistent with what Kazantzidis 
et al. (2004) found. However, we argue that this perturba- 
tion is acceptably small for our purposes. This is because the 
effect on the central dark matter profile is really quite mod- 
erate, and occurs in a radial region that is already dominated 
by baryonic physics when a disk component is included as 
well. In fact, our disk galaxy models show remarkable little 



secular relaxation when started from initial conditions con- 
structed in the above fashion, even for high gas fractions. 
We here benefit from the included treatment of gas pressure 
forces and the improved potential computation compared to 
previous methods (Hernquist, 1993a; Springel, 2000). 

2.4 Galaxy model parameters 

In summary, we specify a disk galaxy model with the fol- 
lowing parameters: 

• The total mass is given in terms of a 'virial velocity'. 
We set Mtot = v^oo/{WGHo). 

• The total mass of the disk is given in terms of a di- 
mensionless fraction rrtd of this mass, i.e. Afdisk = nidAItot- 
Likewise, the mass of the bulge is given in terms of a dimen- 
sionless fraction mb of the total, i.e. Mbuige = mbMtot. The 
rest of the total mass is in the dark matter halo. 

• A fraction /gas of the disk is assumed to be in gaseous 
form, the rest in stellar form. The bulge is taken to be com- 
pletely stellar. 

• The stellar disk is assumed to have an exponential pro- 
file with radially constant vertical scale- height zo- We spec- 
ify the latter in units of the radial disk scale length, and 
typically use zo — 0.1 — 0.2 /i. 

• Similarly, the bulge scale length is specified as a fraction 
of the radial disk scale length. 

• A value for the spin parameter A of the halo is selected. 
Its value effectively determines the radial scale length of the 
disk. 

In Figure 3, we show rotation curves for two galaxy 
models that we use in many of the numerical simulations in 
this study. The two models have a mass comparable to the 
Milky Way, Mtot = 0.98 x lO" fr^MQ, corresponding to a 
virial velocity W200 = 160 kms^^. Their disk mass is specified 
by md — 0.041. The model that includes a bulge (shown on 
the right) has a rotation curve very similar to the one used by 
Hernquist (1993a). Note that the model on the left, with its 
slowly rising rotation curve due to the absence of a bulge (a 
similar model was used, e.g., by Mihos & Hernquist, 1994b), 
is more susceptible to tidal perturbations. 



3 STAR FORMATION AND STELLAR 
FEEDBACK 

Current state of the art simulations lack the dynamic range 
needed to follow the detailed structure of star-forming gas 
on scales characterising entire galaxies. For example, in sim- 
ulations like those presented in the following sections of this 
paper, tens or hundreds of thousands of particles are used to 
characterise the gas in individual disk galaxies, so that each 
particle represents roughly ^ 10^ M© of interstellar mate- 
rial for typical disk gas fractions. However, in order to re- 
solve bound structures in SPH, an object must contain at 
least a number of particles comparable to the number aver- 
aged over when forming smoothed quantities. Here, we set 
this parameter typically to be ~ 64 neighbouring particles, 
so that the simulations cannot resolve bound objects with 
masses smaller than about ^ 6 x 10^ Mq. Simulations with 
total particle numbers 100 times larger than those presented 
here are feasible, but even then we would not able to resolve 
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Figure 4. Effective equation of state for the star-forming gas 
(solid curve) . We show the effective pressure measured in cgs units 
as a function of the gas density in hydrogen atoms per cubic 
centimetre. The vertical dotted line delineates the transition from 
an isothermal gas (dashed) at temperature IC* K to the regime 
governed by an effective pressure in our multiphase model. For 
the particular parameters chosen in this example, the transition 
density lies at 0.128 cm"^. 



self-gravitating structures with masses ^ 10^ M©. Further- 
more, such calculations would be limited to a spatial reso- 
lution of order ^ 20 pc. These mass and spatial scales are 
barely sufficient to resolve individual giant molecular clouds 
and arc inadequate to describe the details of the star- forming 
gas contained within them. 

More seriously, there is at present no complete theory 
for star formation, so the important physics is not well un- 
derstood. Even if the simulations could include the entire 
dynamic range needed to resolve the essentials of star for- 
mation on galactic scales, the description of this process 
would therefore still have to be parameterised using bulk 
averaged quantities. For these reasons, we are motivated to 
adopt an approach in which we discard detailed information 
about star-forming gas in favour of a procedure that cap- 
tures generic consequences of star formation and associated 
feedback effects on resolved scales. In this manner, we do 
not attempt to characterise all aspects of the physics but 
instead seek to identify general trends. 

As an illustration of our procedure, we employ the sub- 
resolution multiphase model for star-forming gas developed 
by Springel & Hernquist (2003a, hereafter SH03). In this 
picture, the ISM is pictured to consist of cold clouds where 
stars form, embedded in a hot, pressure-confining phase. The 
gas is assumed to be thermally unstable to the onset of a 
two-phase medium at densities above a threshold pth- The 
fraction of mass in the two phases at higher densities is set by 
star formation and feedback, evaporation of the cold clouds 
through thermal conduction, and the growth of clouds owing 
to radiative cooling. The rates at which the mass in the two 
phases evolve arc given by eqns. (5) and (6) in SH03. The 
thermal budget of the cold clouds and hot gas is determined 
by the same effects and evolves according to eqns. (8) and 
(9) in SH03. 



Our coarse-grained representation of this model reduces 
largely to two principle ingredients: (1) a star formation pre- 
scription or "law", and (2) an effective equation of state 
(EOS) for the ISM. For the former, we adopt a rate moti- 
vated by observations (Kennicutt, 1989, 1998): 



(22) 



where dp«/di is the star formation rate, /3 depends on the 
initial mass function (IMF) and is the mass fraction of stars 
that die immediately as supernovae, pc is the density of cold 
clouds, and is a characteristic timescale. For a Salpeter 
type IMF we have /3 « 0.1. In our multiphase picture, pc w 
p, where p is the total gas density (sec Fig. 1 in SH03). A 
good match to observations is obtained if is assumed to 
be proportional to the local dynamical time 



U{p) = tl ( ^ ) 



-1/2 



(23) 



where t2 is a constant parameter. 

The effective equation of state (EOS) of the gas in this 
model is given by 



-Pcff = (7^1) {phUh + PcUc), 



(24) 



where pn and Uh are the density and specific thermal energy 
of the hot phase, with similar quantities defined for the cold 
clouds, while 7 = 5/3 is the adiabatic index of the gas. In 
detail, the EOS depends on the values of the parameters in 
our model. 

As described by SH03, the parameter values are con- 
strained by the nature of the cooling curve for the gas and by 
requiring that the EOS be a continuous function of density. 
SH03 also showed that if star formation is rapid compared 
with adiabatic heating or cooling owing to the motion of 
the gas, then the multiphase picture leads to a cycle of self- 
regulated star formation where, in equilibrium, the growth 
of cold clouds is balanced by their evaporation from super- 
nova feedback. Under these conditions and for a given IMF, 
the EOS is determined by three parameters: the constant 
appearing in the star formation law, t1, a normalisation of 
the cloud evaporation rate, Aq, and a supernova "tempera- 
ture", TsN, that reflects the heating rate from a population 
of supernova for the adopted IMF. The requirement that the 
gas be thermally unstable at densities above pth fixes the 
ratio Tsn/^o ~ 10^ while in the self-regulated regime, the 
condition that the EOS be continuous implies that pth de- 
pends mainly on the ratios Ts^/Aq and TsT^/t° (see eqn. 23 
of SH03). 

In Figure 4, we show an example of the form of our effec- 
tive EOS for the parameter choices Tsn = 10*, Aq = 1000, 
and tl = 2.1 Gyrs. For an isolated disk galaxy, these values 
reproduce the observed correlation between star formation 
and gas density (see Figs. 2-3 of SH03). For densities higher 
than Pth, the EOS shown in Figure 4 is fitted to an accuracy 
of 1% by (Robertson et al., 2004) 



logPeff = 0.050 (log nn)^ 
-1- 1.749 logriH 



0.246 (lognn) 
10.6. 



(25) 



The important features of our EOS can be seen from Fig- 
ure 4. If the gas were isothermal at all densities, the ef- 
fective pressure, Peff, would increase linearly with density. 
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as indicated by tlie dashed line for densities higher than 
pth- However, when we account for supernova feedback in 
our multiphase model, the EOS in star-forming gas is stiffer 
(solid curve in Figure 4). Suppose a region at densities lower 
than pth is compressed. If the gas were isothermal, the ef- 
fective pressure would increase proportional to the density. 
However, in our multiphase model, once p > pth, the gas is 
thermally unstable and stars cam form. As the density rises 
further, the star formation rate increases even faster, as does 
the rate of energy injection into the gas from supernova. If 
this feedback energy is retained by the gas, its effective EOS 
of state will be stifFer than isothermal. 

The concept of an "effective equation of state" thus in- 
corporates the impact of local feedback in regulating the 
thermodynamic properties of the ISM. Using this concept, 
we can describe the dynamics of star-forming gas on galactic 
scales and account for the consequences of stellar feedback 
on large scales. For example, the dynamical stability of a 
pure gas disk to axisymmetric perturbations is determined 
by the condition (Toomre, 1964) 



Q 



ttGE 



> 1, 



(26) 



where Cs is the sound speed, k is the epicyclic frequency, 
and E is the gas surface density. This relation represents 
a competition between the stabilising influences of pressure 
(cs) and angular momentum (k) against the destabilising 
effect of gravity (E). If we consider a small region in the 
disk, the appropriate EOS with which to compute Cs is not 
one that captures the detailed structure of the ISM, but an 
effective EOS on the scale of the instability. For parameters 
typical of galactic disks, the fastest growing mode has a 
wavelength that is comparable to the size of the disk (Binney 
& Tremaine, 1987), and so stability is determined by the 
response of the ISM to compressions on this scale. 

In the absence of a complete theory for the ISM, the 
choices for the star formation law and effective equation of 
state are not unique. For example, Barnes (2004) has re- 
cently proposed a star formation prescription that is based 
on the rate at which gas generates entropy in shocks. His nu- 
merical tests indicate that this formalism may yield a bettor 
match to observations of interacting galEixies than a purely 
density-dependent star formation law. 

In our approach, we will consider the macroscopic star 
formation law and effective equation of state as elements 
of our description that can in principle be varied indepen- 
dently to isolate the laxge-scale physics that regulates star 
formation in colliding and merging galaxies. In particular, 
we will also consider simulations where we "soften" the EOS 
by linearly interpolating between that for the standard im- 
plementation of our multiphase model and an isothermal 
EOS. We will denote this softening by a parameter qeos so 
that ^Eos = 1 will correspond to e.g. the "stiff" EOS shown 
by the solid curve in Figure 4, and qeos ~ will correspond 
to the "soft," isothermal EOS, indicated by the dashed curve 
in Figure 4. 

While we do not expect that the effective EOS employed 
here is a precise description of real star-forming gas, the ap- 
proach that wc adopt is flexible and can bo applied to any 
well-specified theory for the ISM. In particular, Springel & 
Hernquist (2003b) have shown that our coarse-graining pro- 
cedure leads to a numerically converged estimate for the 



cosmic star formation history of the Universe (Hernquist & 
Springel, 2003) that agrees well with low redshift observar- 
tions. 

The hybrid multiphase model of SH03 also includes a 
phcnomcnological representation of galactic winds. In this 
picture, some of the energy available from supernova feed- 
back is tapped to drive am outflow if the star formation rate 
exceeds a critical threshold. For simplicity, we do not include 
this additional mode of feedback in the present study. 



4 BLACK HOLE ACCRETION AND AGN 
FEEDBACK 

It is now widely believed that black hole growth and associ- 
ated feedback energy from this process may be important for 
a variety of phenomena related to the evolution of ordinary 
galaxies as well as unusual behaviour in quasars, starbursts, 
and AGN. It has been suggested that the Mbh — cr rela- 
tion may arise if strong outflows arc produced in response 
to a major phase of black hole accretion, which via their 
interaction with the surrounding gas would inhibit any fur- 
ther accretion and hence black hole growth (Silk & Rees, 
1998; Fabian, 1999; King, 2003; Wyithe & Loeb, 2003). In- 
deed, X-ray observations of a number of quasars (mostly 
broad absorption line systems) reveal signiflcant absorption, 
implying large outflows with a kinetic power corresponding 
to a significant fraction of the AGN bolomctric luminosity 
(Chartas ct al., 2003; Crenshaw et al., 2003; Pounds et al., 
2003). In the case of radio- loud QSO, there is also evidence 
that up to half of the total power is injected in the form 
of jets (Rawlings & Saunders, 1991; Tavecchio et al., 2000). 
Inevitably, such outflows must have a strong impact on the 
host galaxies. Models of interacting and merging galaxies 
should, therefore, account for feedback from both star for- 
mation and accretion onto central, supermassive black holes. 
Note that black hole (BH) accretion is expected to be re- 
sponsible for the majority of the BH mass growth and to 
provide sufficient energy supply for driving associated out- 
flows. 

As with star formation, current numerical simulations 
cannot resolve the properties of the accretion flow around 

nuclear black holes on galactic scales. For example, consider 
a black hole of mass Mbh accreting spherically from a sta- 
tionary, uniform distribution of gas whose sound speed at 
inflnity is Coo- The gravitational radius of influence of the 
black hole is then (Bondi, 1952) 

re = —2 — . (27) 

For a black hole of mass Mbh ~ lO'' M© in a gas with a 
sound speed Coo = 30 km/sec, this is numerically 



rB = 50 pc 



Mb 



lO^'M© ) \ 30 km/sec 



(28) 



The Schwarzschild radius of a black hole of this mass is 



2GMbh 



= lO-'^pc 



Mbh 

10^ Mr: 



(29) 



Our understanding of the nature of accretion onto super- 
massive black holes is sufficiently poor that it is not clear 
what range in spatial scales would be required to obtain an 
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accurate description of the impact of black hole growth and 
feedback on galactic scales. However, given the relatively 
poor resolution that can be achieved in simulations like those 
presented here, it is clear that they cannot represent the full 
complexities of this process in any detail. By analogy with 
out treatment of star formation and supernova feedback, we 
are hence led to adopt a coarse-graining procedure. 

In what follows, we use an effective sub-resolution model 
to characterise the growth of supermassive black holes in 
galactic nuclei and the consequences of feedback from ac- 
cretion on surrounding gas. Technically, we represent black 
holes by coUisionless particles that can grow in mass by ac- 
creting gas from their environments. A fraction of the radia- 
tive energy released by the accreted material will be assumed 
to couple thermally to nearby gas and influence its motion 
and thermodynamic state. 

Like our procedure for coarse-graining the ISM, our 
method is flexible and can be applied to any model for 
black hole accretion. As a starting point, we relate the (un- 
resolved) accretion onto the black hole to the large scale (re- 
solved) gas distribution using a Bondi-Hoyle-Lyttleton pa- 
rameterisation (Bondi, 1952; Bondi & Hoyle, 1944; Hoyle & 
Lyttleton, 1939). In this description, the accretion rate onto 
the black hole is given by 



Mb 



(ci-f 1)2)3/2 



(30) 



where p and Cs are the density and sound speed of the gas, 
respectively, a is a dimensionless parameter, and v is the 
velocity of the black hole relative to the gas. We will also 
assume that the accretion is limited to the Eddington rate 



Me 



47r G Mbh rrip 



(31) 



where rrip is the proton mass, ctt is the Thomson cross- 
section, and Er is the radiative efficiency. The latter is related 
to the radiated luminosity, Lr and accretion rate, Mbh, by 



(32) 



Mbhc2 



i.e. it simply gives the mass to energy conversion efficiency 
set by the amount of energy that can be extracted from the 
innermost stable orbit of an accretion disk around a black 
hole. For the rest of this study, we adopt a fixed value of 
Er = 0.1, which is the mean value for radiatively efficient 
Shakura & Sunyaev (1973) accretion onto a Schwarzschild 
black hole. We ignore the possibility of radiatively inefficient 
accretion phases. The accretion rate is then 



Mbh = min(MEdd, Mb) 



(33) 



We will assume that some fraction tt of the radiated 
luminosity Lr can couple thermally (and isotropically) to 
surrounding gas in the form of feedback energy, viz. 



ef Li 



Ef MbYI C 



(34) 



Characteristically we take ef ~ 0.05, so that ~ 0.5% of the 
accreted rest mass energy is available as feedback. This value 
fixes the normalisation of the Mbh — c" relation, and brings it 
into agreement with current observations (Di Matteo et al., 
2004b). 



Figure 5. Evolution of the gaseous disk in an isolated galaxy 
model with a bulge component. A fraction /gas = 0.1 of the disk 
mass is in gas. The remaining 90% is in old stars. The panels show 
a face-on projection of the gas in the disk, and measure 30 h~^kpc 
on a side. The colour-coding indicates both gas density (in terms 
of intensity) and local gas temperature (in terms of colour hue). 
Time in Gyrs is indicated in the upper left corner of each frame. 



5 NUMERICAL APPROACH 

Nothing in our formulation of supernova and black hole 
feedback is tied to a particular numerical scheme for solv- 
ing the dynamical equations for the gas and coUisionless 
matter. Our sub-resolution models and coarse-graining pro- 
cedure could be incorporated into either particle- or grid- 
based methods. In our simulations, we employ an N-body 
algorithm for the coUisionless material which, in our case, 
includes dark matter, stars, and black holes. For the hy- 
drodynamics, we employ a smoothed particle hydrodynam- 
ics (SPH) code which represents fluids elements by parti- 
cles (Lucy, 1977; Gingold & Monaghan, 1977; Monaghan, 
1992). In SPH, fluid properties at a given point are esti- 
mated by local kernel-averaging over neighbouring particles, 
and smoothed versions of the equations of hydrodynamics 
are solved for the evolution of the fluid. 

The particular code we use is an improved and updated 
version of GADGET (Springel et ah, 2001). The code im- 
plements a TreeSPH algorithm (Hernquist & Katz, 1989) 
where gravitational forces are computed with a hierarchi- 
cal tree method (Barnes & Hut, 1986) and SPH is used for 
the hydrodynamics. In the new version GADGET2 employed 
here, the hydrodynamical equations are solved using a fully 
conservative technique (Springel & Hernquist, 2002), which 
maintains strict entropy and energy conservation even when 
smoothing lengths vary adaptively (Hernquist, 1993b). 

Besides gravity and hydrodynamics, the code follows 
radiative cooling processes of an optically thin primordial 
plasma of helium and hydrogen in the presence of an UV 
background (Katz et al., 1996). Star formation and the dy- 
namics of the highly overdense gas are treated with the mul- 
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tiphase model described earlier. Independent star particles 
are stochastically spawned out of the gas phase (Springel & 
Hernquist, 2003a), thereby avoiding artificial coupling be- 
tween gaseous and coUisionless matter. 

Black holes arc represented by coUisionless "sink" parti- 
cles which only feel gravitational forces. For these particles, 
we compute estimates for the gas density in their local envi- 
ronments, in the same fashion as it is done for normal SPH 
particles. Similarly, we also determine the average gas tem- 
perature in the local SPH smoothing environment around 
black hole particles, as well as the gas bulk velocity relative 
to the black holes. These quantities are then used to esti- 
mate the black hole accretion rates, based on the equations 
specified in the previous section. 

To implement the actual accretion, we follow a similar 
stochastic approach as it is applied for regular star forma- 
tion. To this end, we compute for each gas particle j around 
a black hole a probability 



Pi = 



WiMBnAt 



(35) 



for being absorbed by the BH. Here Mbh is the BH accre- 
tion rate. At is the timestep, p is the gas density estimated 
at the position of the black hole, and wj is the kernel weight 
of the gas particle relative to the BH. We then draw ran- 
dom numbers Xj uniformly distributed in the interval [0, 1[ 
and compare them with the Pj. For Xj < pj, the gas parti- 
cle is absorbed by the black hole, including its momentum. 
On average, this procedure ensures that the BH particle ac- 
cretes the right amount of gas consistent with the estimated 
accretion rate Mbh- 

However, because the accretion rate depends sensitively 
on Mbh, this procedure would be quite noisy under condi- 
tions of poor resolution, where accreting a single gas particle 
can change the BH particle mass by a substantial factor. We 
circumvent this problem by giving the sink particle an in- 
ternal degree of freedom that describes the BH mass in a 
smooth fashion. The value of this variable represents the 
BH mass for the case of ideal resolution where the gas mass 
is not discretised. The real dynamical mass of the sink par- 
ticle tracks this internal mass closely, albeit with stochastic 
fluctuations around it. For this reason, we use the internal 
mass for computing the accretion rates, while the growth of 
the hole is followed both in terms of the internal mass and 
the actual dynamical mass. The latter increases in discrete 
steps when whole gas particles are absorbed, but follows 
the internal variable in the mean, with fluctuations that be- 
come progressively smaller for better resolution. Using this 
procedure we can reliably follow the growth of black holes 
in terms of their internal mass variable even in halos with 
just a few hundred particles, which is particularly important 
for cosmological simulations, where the earliest generations 
of galaxies are typically not very well resolved. Note that 
we conserve momentum when gas particles are absorbed by 
BH sink particles. If the BH is moving relative to the gas 
and has a high accretion rate, this can effectively act like a 
friction force which reduces the relative motion. 

Together with the accretion, we compute the rate of 
feedback associated with the black hole growth. This energy 
is added kernel-weighted to the thermal reservoir of the gas 
in the local environment around the black hole. Note that if 
the local cooling rate of the gas is not high enough to radiate 



away all of this energy on a short timescale, an increase 
of the gas sound speed occurs which can then throttle the 
accretion in the Bondi-dominated regime. 

In the final stages of a galaxy merger, the cores of the 
galaxies coalesce to form a single dark matter halo, and a 
single stellar system. Presumably, this also means that a cen- 
tral binary system of two supermassive black holes is formed, 
which may subsequently harden and eventually lead to phys- 
ical collision and merger of the black holes themselves. How- 
ever, it is a controversial question how long it would take to 
harden the black hole binary by stellar-dynamical (Makino 
& Funato, 2004) or hydrodynamical processes (Escala et al., 
2004) until finally gravitational wave emission becomes im- 
portant, leading to quick orbital decay. If black hole merger 
processes were inefficient, it should be a common situation 
that a third supermassive black hole is brought in by the 
next galaxy merger. This black hole could then interact with 
the binary and lead to sling-shot ejection of one of the holes, 
with the remaining pair being ejected in the opposite direc- 
tion. As this may seriously impair the growth of black holes 
to the large masses that axe observed, the existence of super- 
massive black holes may be taken as circumstantial evidence 
that BH binaries probably do merge. We therefore assume 
that two black hole particles merge instantly if they come 
within each others's smoothing lengths, and their relative 
velocities are at the same time smaller than the local sound 
speed. Note that as with the accretion itself, we lack the 
dynamic range in simulations of whole galaxies to study the 
hardening process of black hole binaries directly. 



6 ISOLATED DISK GALAXIES 

We have run a series of simulations of isolated galaxies us- 
ing the methods described in Section 2 to initialise them. 
In the following, we examine their stability using models 
with and without bulges, having the rotation curves shown 
in Figure 3. The structure of these models was motivated 
by properties of gala^xies like the Milky Way and to enable 
us to compare our results with earlier work, such as Hern- 
quist & Mihos (1995) and Mihos & Hernquist (1996). Here, 
we focus on the effect of varying the stiffness of the EOS 
and the gas fraction in the disk. A significant advantage of 
our formalism is that we are now able to construct stable 
equilibrium models with a larger gas fraction than possible 
in previous studies. 



6.1 Stability of models without black holes 

In order to facilitate comparison between different cases, 
in this section we fix the parameters associated with our 
multiphase model and vary only the disk gas fraction and 
the stiffness of the effective EOS. For definiteness, we take 
1° = 8.4 Gyrs, Ao = 4000, and Tsn = 4 x 10*. This choice for 
t" is a factor of four larger than that employed by SH03 to 
match the Kennicutt Law. However, we find that for our 
isolated galaxies, t'i — 2.1 Gyrs yields global star forma- 
tion rates ^ 4M0/yr for gas fractions of 10% and struc- 
tural properties similar to the Milky Way. In our Galaxy, 
the inferred average star formation rate is only 1 M© /yr. 
Increasing t° by a factor of four gives therefore better agree- 
ment with the long gas consumption timescale inferred for 
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Figure 6. Stability analysis of disk galaxies as a function of gas fraction in the disk, /gas, and stiffness of the effective EOS. Each 
panel shows the face-on distribution of gas obtained at time t = 0.3 Gyrs after the start of the simulations, roughly corresponding to 
one rotation period at three exponential scale-lengths. The gas fraction is constant along columns, while each row shows the results for 
a fixed EOS. In all cases, the initial galaxy models had the same structure, except for a different partitioning of the mass in the disk 
between a gaseous and a stellar component. The galaxies in these examples did not include stellar bulges, but the results for galaxy 
models with bulges are qualitatively very similar. 



the Galaxy. This also simplifies a comparison with Hern- 
quist & Mihos (1995) and Mihos & Hernquist (1996), who 
chose parameters in their star formation prescription to give 
global star formation rates ~ 1 Mq /yr for galaxies like the 
Milky Way with 10% of their disk mass in gas. 

Once t'^ is fixed, we select Aq and Tsn as above so that 
the critical density for the gas to be thermally unstable, pth, 
implies a critical projected gas surface density for star for- 
mation to occur that is similar to observations (Kennicutt, 
1989, 1998). Since pth in our formalism depends mainly on 
the ratios Tsn /Aq and Tsn /t^l , a factor of four lengthening of 
relative to SH03 requires a factor of four increase in Tsn 
and Ao for a critical gas surface density ~ IOMq/pc"'^ (see 
Figs. 2-3 of SH03). As we indicated earlier, the effective EOS 
depends primarily on ratios between these three parameters. 



Scaling t", Tsn, and by the same factor implies that the 
EOS for our isolated disks is the same as the "stiff" case in 
Figure 4, only the gas consumption timescale is larger. In 
the stability analysis shown below, we hold these parame- 
ters fixed, so that the threshold density for star formation is 
the same, but we artificially vary the stiffness of the EOS by 
linearly interpolating between the stiff and isothermal limits 
in Figure 4. EOS softenings of qeos = 1 and qeos ~ refer 
to the stiff and isothermal cases, respectively. We also vary 
the fraction /gas of the disk mass in gas, but keep the total 
disk mass constant. 

As a typical example. Figure 5 shows the time evolution 
of the projected gas density in an isolated galaxy model with 
a bulge. The gas fraction is /gas = 0.1 and the EOS soften- 
ing parameter is qeos = 0.5. Initially, the particles randomly 
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Figure 7. Evolution of the global star formation rate in galaxies 
without bulges, as a function of the gas fraction and EOS used. 
The star formation rate in solar masses per year is scaled by the 
disk gas fraction raised to the 1.4 power (as in the Kennicutt 
Law). Each panel shows the evolution for the six values of /gas 
shown in Figure 6, with lighter colours corresponding to larger gas 
fractions. From top to bottom, the panels differ in the employed 
EOS for the star forming gas. Results are given for the stiff EOS of 
the multiphase model (<?eos = !)> an intermediate case between 
this model and an isothermal EOS (^eos = 0.25), and a 'soft' 
model with a nearly isothermal EOS (gEOS = 0.05). 

sample the exponential surface density distribution, so there 
are no coherent structures present in the disk. However, a 
steady spiral pattern develops within a fraction of a rotation 
period owing to the action of swing amplification (Toomre, 
1981). This behaviour is similar to that reported in previ- 
ous numerical studies of disk structure. For the parameter 
choices in this example, the disk is stable for many rotation 
periods and does not evolve strongly. The gas is converted 
into stars at a relatively low rate according to our star for- 
mation prescription but, otherwise, the disk, bulge, and halo 
do not evolve. 

However, strong, unstable evolution is possible if the 
gas fraction is large and the EOS is too soft. In Figure 6, we 
compare the distribution of gas in model galaxies without 
bulges after approximately one rotation period of evolution. 
The gas fraction is increased from /gas ~ 0.1 to /gas =0.99 



Figure 8. Same as in Figure 7, but for galaxy models that include 
bulges. 

along rows (left to right) while the EOS is softened from 
9EOS ~ 1 (top) to Qeos = 0.05 (bottom) along columns. 
Some interesting features and trends are apparent from Fig- 
ure 6. 

First, for a relatively stiff EOS, gEOS ^ 0.3, increasing 
the gas fraction actually suppresses the amplitude of spiral 
structure in the disk. This is because in our full multiphase 
model, the effective temperature of the gas is Tofi ^ 10^ K 
over most of the disk for q — 1, and so the gas is dynam- 
ically hotter than the old stars. Increasing the gas fraction 
in this case makes the disk more stable and less susceptible 
to amplifying non-axisymmetric patterns. Thus, in the top 
row of Figure 6, the spiral structure which is apparent for 
/gas = 0.1 (top left) is mostly washed out for /gas = 0.99 
(top right). This tendency does not occur if the EOS is soft 
because then the effective temperature is such that the gas 
is dynamically colder than the old stars. For ^eos < 0.125 
(bottom two rows), this results in unstable behaviour if the 
gas fraction is large. 

Second, the set of models shown in Figure 6 clearly de- 
lineates stable cases from unstable ones. For our full, stiff 
EOS (qeos = 1) even pure gas disks are stable (top row). 
This is also true for intermediate cases with qeos ~ 0.5. 
However, if the effective temperature of the gas is lower than 



© 0000 RAS, MNRAS 000, 000-000 



Modeling feedback from stars and black holes in galaxy mergers 13 



that of the stars, very gas-rich disks are unstable and frag- 
ment within a rotation period. Thus, for example, with a 
nearly isothermal EOS (qeos ~ 0.05, bottom row), disks 
with a gas fraction /gas ^ 0.3 are already unstable (bottom 
row). 

The previous studies by Hernquist & Mihos (1995) and 
Mihos & Hernquist (1996) treated the ISM as an isothermal 
gas and included only weak feedback in the form of a small 
input of random kinetic energy. These authors were unable 
to construct stable models with large gas fractions owing to 
disk fragmentation like that shown in Figure 6. Supernova 
feedback in our multiphase model pressurises the gas making 
the effective EOS stiffer, stabilising the disk. 

A similar separation between stable and unstable be- 
haviour is seen in our galaxy models that include bulges. We 
illustrate this in another manner in Figures 7 and 8 where 
we show the evolution of the global star formation rates in 
the disks for models without and with bulges, respectively. 

The behaviour in Figures 7 and 8 clearly separates sta- 
ble and unstable models. In unstable disks, fragmentation 
leads to a runaway growth in density perturbations, yield- 
ing unsteady evolution in the global star formation rate. 
Thus, all our galaxies are stable for stiff equations of state 
with gEOS > 0.5. Models with a bulge are stable for all gas 
fractions with gEOS ~ 0.25, however those without a bulge 
(Figure 7) are unstable with this EOS for /gas ^ 0.6. Softer 
equations of state yield unstable behaviour for galaxies with- 
out bulges for /gas ~ 0.4. For models with bulges, instability 
sets in at /gas > 0.8 for qeos = 0.125 and /gas ~ 0.6 for 
9EOS = 0.05. The scalings with gEOS and /gas are broadly 
consistent with the Toomre criterion. 



6.2 Models with black holes 

We have also run models similar to those described in Sec- 
tion 6.1, but where we added a sink particle to represent 
a black hole and allowed it to accrete gas from the galaxy, 
increasing its mass and adding a source of feedback energy 
in addition to supernovae. The examples we discuss in this 
section are not intended to provide a physical picture for 
the growth of supermassive black holes in isolated galaxies 
because we have made no attempt to tie the origin of the 
black holes to galaxy formation. Nevertheless, we can use 
these simulations to identify some generic features of our 
description of the interaction between black holes and the 
ISM that are relevant for dynamically evolving situations, 
like galaxy mergers. 

The coupling between black hole accretion and sur- 
rounding gas has several implications. First, heating from 
AGN feedback energy can drive outflows from galaxies if 
the accretion rate is sufficiently high. This is shown in Fig- 
ure 9, where the effect of black hole feedback is illustrated 
for an isolated galaxy with a bulge. For the low accretion 
rates relevant to the times indicated, the deposition of ther- 
mal energy into the gas in the vicinity of the black hole can 
drive a weak wind perpendicular to the plane of the disk, 
as indicated by the velocity vectors in this figure. These 
outflows are reminiscent of those in our supernova-driven 
winds (Figs. 5-8 in SH03). More powerful outflows can be 
produced during mergers when gas is driven into the centre 
of the remnant, fueling strong nuclear accretion. 

As implied in Section 4, our description of black hole ac- 
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Figure 10. Evolution of the black hole mass in simulations of an 
isolated gas-rich galaxy. We show the black hole mass with (solid) 
and without (dashed) thermal feedback, respectively. 



cretion yields several distinct phases of evolution, depending 
on the properties of the gas near the black hole. For accre- 
tion rates lower than Eddington, the black hole mass will 
evolve according to 



MBH{t) = 



Mo 



(36) 



(37) 



I - xMot' 

where Mo is the initial black hole mass and 

_ 47r aG^ p 
^ ~ (c? + t;2)3/2 ■ 

This result is valid only when the properties of the gas near 
the black hole and the relative velocity do not evolve signif- 
icantly with time. When the accretion rate is higher than 
Eddington, the black hole grows exponentially with time. 



Msnit) =A/bh(0) exp(^) , 
where the Salpeter time. 



ts 



47r G m„ 



(38) 



(39) 



depends only on physical constants and on the radiative 
efficiency. For a radiative efficiency of 10%, ts ~ 4.5 x 10^ 
years. 

In the absence of feedback effects, a low-mass black hole 
would grow first according to equation (36). Once the black 
hole becomes sufficiently massive, given the properties of 
the surrounding gas, it then enters a period of exponen- 
tial growth, described by equation (38). Note that the first, 
Bondi-limited growth can be very slow if the initial black 
hole seed mass is small, or if a small value for the coefficient 
a is selected. We typically start our calculations with black 
hole seed masses of order 10^ Mq and set a large enough to 
allow a hole in a gas-rich environment to reach the Edding- 
ton regime in about ~ 0.5 Gyr. 

We illustrate this behaviour in Figure 10, which shows 
the growth of black holes in isolated, gas-rich galaxies, with 
and without the impact of accretion feedback on the gas. 
If feedback is neglected, the black hole grows in a manner 
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Figure 9. Edge-on view of the gas in an isolated disk galaxy with a bulge that hosts a growing, supermassive black hole in its nucleus. 
The two panels of side-length 32 h~^kpc show snapshots at times 0.7 Gyrs (left) and 1.4 Gyrs (right) after the beginning of the simulation. 
At the low accretion rates present in this quiescent galaxy, feedback energy from the black hole can nevertheless drive a weak wind into 
the halo. 



consistent with the expectations noted above, first according 
to equation (36) and later exponentially. When the effects of 
feedback are included, however, the growth is self-regulated 
and leads to a saturated final state, that depends on the 
initial gas content of the galaxy. The self-regulation occurs 
because as the accretion rate increases, the amount of energy 
available to affect nearby gas also grows. If this material is 
then heated and dispersed, the supply of gas to grow the 
black hole will decline as accretion is shut off. 



7 MAJOR MERGERS 

The approach we have adopted for handling various forms 
of feedback makes it possible to explore physical effects that 
were not accessible previously. For example, Hernquist & 
Mihos (1995) and Mihos & Hernquist (1996) were able to 
model galaxies with only small gas fractions because their 
treatment of feedback was insufficient to stabilise highly gas- 
rich disks. With our method, we can examine mergers be- 
tween galaxies with a very large gas content, as might be 
relevant to objects at high redshifts. Also, by including pro- 
cesses related to black hole growth and accretion, we can 
study the interplay between stellar and AGN feedback that 
is likely to be important for understanding the nature of 
ultraluminous infrared galaxies (ULIRGs). 

In what follows, we illustrate some of the consequences 
of our models for stellar and AGN feedback by looking at 
major mergers of disk galaxies. A representative case is given 
in Figure 11, which shows the distribution of gas during the 
merger of two equal-mass disk galaxies from a prograde, 
parabolic orbit. As the disks pass by one another for the 
first time {t = 0.2 — 0.3), strong tidal forces yield extended 
tails and bridges. Gas is shocked at the interface between 
the two galaxies and the tidal response drives gas into the 



central regions of the two galaxies. When the galaxies finally 
coalesce {t = 1.2 — 1.5), much of the remaining gas is either 
converted into stars during an intense burst, or shock-heated 
to temperatures characteristic of the virial temperature of 
the remnant. 



7.1 Mergers without black holes 

7.1.1 Comparison with previous work 

One of our goals is to see to what extent the strength of 
stellar feedback influences starbursts during mergers. As a 
starting point, we compare our results to those obtained in 
earlier studies, in particular, to the work of Hernquist & 
Mihos (1995) and Mihos & Hernquist (1996). These authors 
employed a highly simplified treatment of star formation by 
modeling the ISM as an isothermal gas with a temperature 
TiSM ~ 10^ K, and including a weak form of supernova feed- 
back by injecting small amounts of kinetic energy into the 
gas. The amplitude of this feedback was adjusted so that iso- 
lated disks forming stars quiescently had gas profiles similar 
to those observed. 

We have modified our simulation code to mimic the 
star formation and feedback algorithm of Mihos & Hern- 
quist (1994b) in an attempt to reproduce their results. In 
Figure 12, we show the evolution of the star formation rate 
in prograde mergers of equal-mass disk galaxies with and 
without bulges. We have selected the parameters and galaxy 
models to match the results of Mihos & Hernquist (1996) as 
closely as possible. The results in Figure 12 can be compared 
directly to e.g. Fig. 5a in Mihos & Hernquist (1996). 

The morphologies of the merging galaxies roughly fol- 
lows the time sequence in Figure 11, with only a mod- 
est dependence on their structural properties. When the 
galaxies first pass by one another, relatively strong star- 
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Figure 11. Evolution of tlie gas distribution in a major merger of two disk galaxies. Colour hue encodes gas temperature, while brightness 
indicates gas density. Each frame measures 50/i~^kpc on a side, and the corresponding time of each image is given by the labels. In this 
simulation, a slightly softer EOS (gEOS = 0.25) than in our default multiphase model was used. 



bursts are triggered in each disk in the case when bulges are 
not present. Galaxies with bulges experience much weaker 
starbursts during first passage. Mihos & Hernquist (1996) 
showed that this difference arises from the relative ease with 
which the disks can amplify m = 2 disturbances in response 
to tidal forcing. During final coalescence, a much strong star- 
burst is excited in the models with bulges because much of 
the gas in the bulgeless galaxies has already been consumed 
by this time. These trends agree well with those identified 
by Mihos & Hernquist (1996). 

The amplitudes of the starbursts shown in Figure 12 are 



comparable to, but somewhat stronger than those shown 
in Fig. 5a of Mihos & Hernquist (1996), and considerably 
stronger than those found by Barnes (2004) (see e.g. his 
Fig. 5). By varying the parameters of our simulations, we 
have found that the evolution of the star formation rate 
in mergers with only relatively weak feedback is sensitive 
to numerical resolution. If the star formation rate is tied to 
gas density, the amplitudes of merger-induced starbursts de- 
pend on the compressibility of the gas, which is infiuenced 
by both the stiffness of the EOS, as well as dynamic range in 
resolution of the numerical algorithm. We have verified that 
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Figure 12. Evolution of the star formation rate in major mergers 
of disk galaxies, using a description of the ISM that mimics the 
approach of Hernquist & Mihos (1995) and Mihos & Hernquist 
(1996). We show results for models without (left) and with (right) 
bulges. The galaxies have the rotation curves shown in Figure 3, 
and the disks initially consisted of 10% gas. 



diflterences in resolution are mainly responsible for the resid- 
ual discrepancies between the star formation rates shown in 
Figure 12 and those reported by Mihos & Hernquist (1996) 
and Barnes (2004). 



7.1.2 Multiphase treatment of ISM 

As a next step, we have performed a grid of merger simu- 
lations using the multiphase treatment of supernova feed- 
back described in Section 3. We have varied the structural 
properties of the galaxies, the orbit of the mergers, and the 
parameters in our multiphase model. Below, we limit the 
discussion mainly to the consequences of variations in the 
fraction of gas in each galaxy and the strength of feedback 
to highlight the new types of outcomes that are possible. 

In Figure 13, we show the evolution of the star for- 
mation rate from one of our simulations, a prograde, major 
merger of two disk galaxies without bulges. In this case, each 
disk initially consisted of /gas ~ 99.9% gas and only 0.1% 
of old stars. The EOS was relatively stiff with a softening 
parameter of gEOS = 0.5. (That is, the EOS was linearly 
interpolated to lie midway between the curves in Figure 4 
for p > pth.) Even with this extreme gas content, the model 
disks are stable in isolation, owing to the pressurisation of 
the gas from supernova feedback. The parameters of this 
simulation are deliberately chosen to be extreme to simplify 
the discussion, but qualitatively similar results are obtained 
in more general situations. 

In the example shown in Figure 13 the star formation 
rate peaks at roughly 500 M0yr~'^, considerably higher than 
in earlier simulations of equal-mass mergers of disks, owing 
to the large gas content of each galaxy. Star formation rates 
at this level are similar to those inferred for systems at high 
redshift such as Lyman-break galaxies and SCUBA sources. 



Figure 13. Evolution of the star formation rate in a major merger 
of two disk galaxies without bulges, and with disks that initially 
consisted of 99.9% gas. 

suggesting that some of these objects may result from merg- 
ers of gas-rich disks. 

Figure 13 also demonstrates that while the general con- 
clusions of Mihos & Hernquist (1996) are obtained under 
broader circumstances, some of their results clearly depend 
on the strength of feedback from star formation. In mod- 
els with weak supernova feedback, as in Figure 12, galaxies 
without bulges experience stronger starbursts at first pas- 
sage than during final coalescence. The evolution shown in 
Figure 13 is rather difi'erent, with the strongest starburst ac- 
companying the final merger. While the merging galaxies in 
this example still develop a bar-like structure in response to 
tidal forcing, the gas is sufficiently hot dynamically that it is 
not strongly compressed during this phase of evolution, lim- 
iting the strength of the initial starburst. The more intense 
starburst during the final merger results from strong shock 
compression of the gas and gravitational torquing when the 
galaxies coalesce, in a manner reminiscent of the mergers 
with bulges examined by Mihos & Hernquist (1996); Barnes 
& Hernquist (1991, 1996). 

An interesting aspect of the simulation used in Figure 13 
concerns the structure of the merger remnant. A consider- 
able amount of gas is left behind following the merger. Much 
of it settles into an extended disk, surrounding a relatively 
compact, rapidly rotating remnant of newly formed stars. 
Overall, the remnant is morphologically closer to a spiral 
galaxy with a bulge than to an elliptical. Further implica- 
tions and a detailed analysis of this simulation are presented 
in Springel & Hernquist (2004). 

7.2 Mergers with black holes 

An even greater diversity of outcomes results during a 
merger if the impact of central, supermassive black holes is 
included in addition to feedback from star formation. As we 
discuss in Di Matteo et al. (2004b) and Springel et al. (2004) , 
mergers of galaxies that involve star formation, supernova 
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feedback, black hole growth and accretion, and AGN feed- 
back yield a generic time history in which the tidal inter- 
action drives gas into the centre of the remnant, fueling a 
starburst and triggering the rapid growth of the central black 
hole(s). Feedback from accretion onto the black hole(s) in- 
fluences the thermodynamic properties of the surrounding 
gas, eventually expelling it from the inner regions of the 
remnant, terminating the episode of rapid star formation 
and self-limiting the growth of the black hole(s). The star- 
burst and AGN activity are coeval, but offset in time owing 
to the detailed form of the response of the gas to the feed- 
back. Thus, while we expect starbursts and AGN activity 
to be correlated in this picture, the remnant will be primar- 
ily seen in different stages as the merger progresses and will 
evolve from one type of object into the other, depending on 
the observed wavelength regime. 

In Figure 14 we show three snapshots from a simulation 
that includes our model for the growth of black holes. In this 
example, the galaxies included bulges and the disks were 
10% gas. Following the first passage (shown on the left), but 
before the galaxies coalesce, the disks are distorted by their 
mutual tidal interaction, but only a relatively weak starburst 
is triggered, because the bulges stabilise the disks against 
a strong m = 2 response. The black holes do not accrete 
significantly during this phase of evolution, and so neither 
galaxy would be observable as an AGN. When the galaxies 
merge but before the gas has been consumed (middle), tidal 
forces trigger a nuclear starburst and fuel rapid growth of the 
black holes. The peak in the starburst occurs near the time 
indicated in Figure 14, but the peak in the black hole growth 
is offset, owing to the delayed action of AGN feedback on the 
gas. In principle, the system should now be both a starburst 
and an AGN, but it is likely that the quasar-activity would 
be obscured by surrounding gas and dust, at least during 
the beginning of the burst. Only during the final stages of 
evolution could the remnant be visible as an AGN, when 
outflows remove the dense layers of gas around it. 

As the remnant settles into a relaxed state, star forma- 
tion and black hole accretion are rapidly quenched by the 
expulsion of gas from the centre, as a consequence of AGN 
feedback. The remnant would no longer be observable as ei- 
ther a starburst or an AGN and will age rapidly, quickly 
resembling an evolved red galaxy (Springel et al., 2004). At 
this point the black hole growth ceases and the black hole 
mass attained is consistent with those expected from the 
observed Mbh — o" relation (Di Matteo et al., 2004b). 

The interplay between black hole growth and the 
physics of the ISM in simulations like that shown in Fig- 
ure 14 has significant implications for the observed proper- 
ties of these systems. In Figure 15, we show the evolution 
of the intrinsic star formation rates in simulations of galaxy 
mergers with and without bulges and including or neglect- 
ing central black holes. This figure shows that the presence 
of the black holes and the feedback energy derived from 
their growth can dramatically alter the nature of the star- 
burst resulting from a merger. The peak amplitude of the 
starburst in the merger which included bulges is lowered by 
about a factor of five, owing to AGN feedback. We note, 
however, that the observed emission from such an object 
would include both the radiation from the starburst as well 
as the energy output from black hole accretion. If most of 
the AGN emission is reprocessed into optical or infrared ra- 
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Figure 15. Star formation rates in mergers of galaxies without 
(left) and with (right) bulges. In each panel, the two curves show 
the outcome with and without black holes. 

diation by surrounding gas and dust before it escapes from 
the inner regions of the remnant, it could be incorrectly at- 
tributed to the starburst, rather than being associated with 
a buried quasar. A more complete understanding of the im- 
plications of these results will require a detailed modeling 
of the spectral energy distribution from the combination of 
the starburst and the obscured AGN. 

We have varied the orbital geometry and structural 
properties of the merging galaxies to see if the results 
described above are generic. While the detailed outcome 
does depend on e.g. the orbit, many of the overall prop- 
erties are relatively unaffected. An example is shown in Fig- 
ure 16 where the evolution of the total black hole mass is 
given for six different orbital configurations. The top three 
panels show prograde-prograde, prograde-retrograde, and 
retrograde-retrograde mergers (left to right), while the bot- 
tom three panels show cases where the disk orientations were 
chosen randomly. While the detailed response varies from 
model to model, the final total black hole mass is insensitive 
to the orbital parameters. In all cases, the final black hole 
mass is ~ 2 — 3 X 10^ M© . This value is consistent with the 
observed correlation between black hole and spheroid mass 
for the remnants produced in these mergers, as shown in 
more detail in Di Matteo et al. (2004b). The simulation re- 
sults are in accord with the view that black hole growth is 
regulated by AGN feedback and the dynamical response of 
gas to this supply of energy. 

We note that in our model the AGN feedback energy 
is deposited isotropically in the gas around the black hole. 
While this appears to be a natural approximation given that 
we cannot resolve the detailed accretion physics close to the 
hole, observations show that black hole outflows are typi- 
cally bi-polar (although not necessarily strongly coUimated) 
and often result from powerful jets. It is possible that our 
isotropic coupling enhances the effects of the feedback, but 
this effect should be compensated in part by our compara- 
tively conservative assumption for ef. We also caution that 
the physical nature of the coupling of AGN feedback to the 
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Figure 14. Merger of two disk galaxies, including the effects of black hole growth and AGN feedback. The images show the gas 
distribution in the two disks at three different times, where colour hue encodes temperature while brightness measures gas density. The 
bottom panels show the time evolution of the accretion rate onto the black holes (top) and the star formation rate (bottom). The red 
symbols in these panels mark the three times shown in the images on top. The first snapshot shows the system just after the first passage 
of the two disks. The second snapshot captures the system when the galaxies are coalescing, at which point the star formation and 
accretion rates peak. Finally, the third snapshot shows the system after the galaxies have fully merged, and most of the mass has settled 
into a slowly evolving remnant. 



surrounding gas is not fully understood yet. In particular, 
there may be important feedback meclianisms other than 
thermal heating. For example, radiation pressure may play 
a very important role as well (Murray et al., 2004). 

The presence of central black holes can also have a sig- 
nificant impact on the stellar density structure of merger 
remnants. Mihos &i Hernquist (1994a) found that their rem- 
nants typically had centrally peaked luminosity profiles from 
the compact stellar cores left behind by starbursts. In our 
new simulations, the remnants have smoother luminosity 
profiles into their centres, in much better agreement with 
observations of elliptical galaxies (Springel, Di Matteo & 
Hernquist, 2004, in preparation). This outcome is driven by 
a combination of the consumption of dense gas by the black 
holes, the dynamical heating of the gas as the black holes 
spiral together and merge, and the dispersal of gas by AGN 
feedback. Our results thus indicate that black hole growth 
may be an integral part of spheroid formation. 



8 CONCLUSIONS 

In this paper, we have introduced a new methodology 
for simultaneously modeling star formation and black hole 
accretion in hydrodynamical simulations of isolated and 
merging galaxies. Our approach uses the concept of form- 
ing "macroscopic", coarse-grained representations of sub- 
resolution physics, allowing us to study the effects of these 
processes on larger, resolved scales. To illustrate the princi- 
pal effects of our methods and of the physics we describe, we 
have focused on simulations of isolated and merging galax- 
ies. Throughout, we have concentrated on the methodologi- 
cal aspects of our work, rather than on an in-depth analysis 
of the physical results of our simulations. We provide the lat- 
ter in related papers (e.g. Di Matteo et al., 2004b; Springel 
et al., 2004). 

For our simulations, we have constructed compound 
galaxy models that consist of a dark matter halo, a rota- 
tionally supported disk of gas and stars, and a central bulge. 
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Figure 16. Evolution of the total black hole mass in merger simulations that differ in their encounter geometries. The top three panels 
show results where the disks are oriented in the orbital plane, while the bottom three panels are for simulations with random disk 
orientations. The angles (8, </>) specify the orientation of the spin vector of one of the disk galaxies relative to the orbital plane (see Due 
et al., 2000, for a sketch of the orbital configuration). 



with independent parameters describing each of the struc- 
tural components. Our approach improves upon previous 
methods by accounting for gas pressure forces, as weU as 
featuring a flexible numerical approach for computing the 
gravitational potential for non-trivial, radially varying mass 
distributions in the disk. As a result, we are able to construct 
very stable disk galaxy models even for high gas fractions. 

In simulations of isolated galaxies, we have emphasised 
the importance of the treatment of feedback processes in 
the dense star-forming ISM for the stability of disk galax- 
ies. Simulations that use a simple isothermal equation of 
state, or a weak form of kinetic feedback, have gas disks 
that are very susceptible to axisymmetric perturbations. As 
a result, galaxies with large gas surface mass density cannot 
be evolved in a stable fashion for a long time in such models. 
Our multiphase model for the ISM encapsulates the effects 
of local supernova feedback in the form of an EOS, which 
is comparatively stiff. When this description for the ISM is 
chosen, the dense gas is stabilised by supernovae pressurising 
the gas, so that much larger gas fractions than in previous 
studies become possible. 

We used simulations of isolated galaxy models to ex- 
amine the properties of the various phases of accretion pos- 
sible in our model for the growth of black holes. In par- 
ticular, we showed how the black hole growth of a small 
seed black hole accelerates in the Bondi regime and even- 
tually reaches Eddington-limited, exponential growth. How- 
ever, the growth process can be slowed down and regulated 



by the feedback energy associated with the accretion. We 
assumed a simple thermal coupling of a small fraction of the 
hole's bolometric luminosity to the surrounding gas; this can 
increase the local sound speed of the gas and thereby reduce 
the Bondi accretion rate, or in extreme cases of high accre- 
tion, generate a powerful, pressure- driven quasar wind. 

Using a number of simulations of major galaxy mergers, 
both with and without black holes, we have explored a few 
basic consequences of our methods. When we model the star- 
forming gas with an isothermal equation of state and weak 
feedback, our results are in good agreement with previous 
work by Mihos & Hernquist (1996). However, our multiphase 
treatment of the ISM allows us to investigate galaxy mergers 
with much higher gas fractions than possible before. Here, 
new types of outcomes are possible for major mergers. For 
example, very high star formation rates of several hundred 
solar masses per year can be reached in gas-rich mergers. 
In merger simulations with black holes we have shown that 
the presence of accreting black holes can dramatically alter 
the dynamics of the merger. The feedback energy associated 
with the growth of the hole owing to the tidally triggered 
inflow of gas has a direct effect on the strength of the nu- 
clear starburst. During the final galaxy coalescence, the BH 
can expel the gas from the centre in a powerful outflow, 
quenching the starburst on a short timescale. The resulting 
elliptical remnant is then gas-poor and reddens quickly. 

Already, these exploratory results show clearly that 
adding an accreting supermassive black hole in the cen- 



© 0000 RAS, MNRAS 000, 000-000 



20 V. Springel, T. Di Matteo, and L. Hernquist 



tre can have a decisive influence on galaxies. In fact, we 
think this indicates that the build-up of supermassive black 
holes and massive galaxies are an intertwined process which 
requires a self-consistent treatment to be meaningfully ad- 
dressed. Our new numerical methodology represents the first 
attempt to do carry out this demanding task in cosmologi- 
es! hydrodynamical simulations, opening up the possibility 
to jointly study the formation of galaxies and the build-up 
of supermassive black holes hosted by them. Our methods 
thus provide a powerful tool to gain new insights into the 
process of hierarchical galaxy formation. 
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